% GAUTIER LE BIHAN - 2020
% Replication files for "Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector 
% Menu-Cost Model ?" Review of Economics and Statistics
%
% Estimation of Aggregate 2-sector CP model (energy + core)
% inputs for IRF Table 4
% 

clear all
close all
clear matrix
clc

addpath('..\..\Utilities')  
load actual_moments_k

load res_estim_prod2nc

p=res_estim_prod2nc
param0=   [p(1)  p(3)  p(5)  ;
           p(2)  p(4) p(6) ]

       
actual_moments=[0.0852  0.7189  0.0489  0.0742 5.5153;
                0.6812  0.6110  0.0314  0.0581  2.8919]
actual_std=[0.0009; 0.0053; 0.0008; 0.0018; 0.2666;
    0.0009; 0.0053; 0.0008; 0.0018; 0.2666];
      %      0.0012 ;0.0021; 0.0001 ;0.0003 ;0.0254];

 actual_var= actual_std.* actual_std*1000;

 for jj=1;

   p01 =param0(1,1);
   p02 =param0(2,1);
   
    mu_c01 =param0(1,2);
    mu_c02 =param0(2,2);
    
    sig_eps_a01=param0(1,3);
    sig_eps_a02=param0(2,3);
    
    %rho_a01 =param0(1,4);
     % rho_a02 =param0(2,4);

weight_j=0.05;


% Adjustement of initial search area (beyond 5 percent)

power_p0 = 1.0;
power_mu_c = 1.;
power_sig_eps_a = 1.;
power_rho_a = 1;

% save SMM

% create vector for initial parameter values
vec0 = [p01; p02; mu_c01; mu_c02; sig_eps_a01; sig_eps_a02]

% Keep track of results
conv_SMM = zeros(0,3);
xopt_SMM = 100000; % start iteration
conv_opt = zeros(3,6);
SMM_count = 0; % counter

% tolerance level
tol_SMM = 5e-3;
scale=[actual_var];


save('SMM','vec0','p01','p02','mu_c01','mu_c02','sig_eps_a01','sig_eps_a02',...
        'power_p0','power_mu_c','power_sig_eps_a','power_rho_a',...
    'conv_SMM','xopt_SMM','conv_opt','SMM_count',...
    'tol_SMM','-append', 'weight_j')
% 'power_phi','phi_init','phi0'

%% Options
%  Display, TolX, TolFun, MaxFunEvals

%% Set calibration target
% std.Y, autocorr(Y), frequency of default:
moment1 = actual_moments(1,1);
moment2 = actual_moments(1,2);
moment3 = actual_moments(1,3);
moment4 = actual_moments(1,4); 
moment5 = actual_moments(1,5);

moment6 = actual_moments(2,1);
moment7 = actual_moments(2,2);
moment8 = actual_moments(2,3);
moment9 = actual_moments(2,4); 
moment10 = actual_moments(2,5);


target = [moment1; moment2; moment3; moment4; moment5; moment6; moment7; moment8; moment9; moment10];
save('SMM','target', 'scale','-append')
%% Start minimization routine 


%options = optimset('MaxIter',30);%, 'TolX', 1e-4,'Tolfun', 1e-3);
options = optimset('TolX', 1e-4,'Tolfun', 1e-3);
options = optimset('MaxIter',20);
options = optimset('MaxFunEvals',100);

[x,fval, exiflag] = fminsearch(@SMMmain,vec0, options);
res_estim_prod2nc(jj,:)=[x' fval exiflag target' ]

 
  save res_estim_prod2nc res_estim_prod2nc ;

end